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ABSTRACT 

In  part  I,  several  simple  Inltlal-value  problems 
based  on  an  equation  of  the  form 

where  ^  is  a  given  function,  are  considered,  and  the 
stability  of  certain  difference  schemes  for  these 
problems  is  studied  numerically.   The  leap-frog  scheme 
(with  and  without  a  dissipative  term),  the  Lax-Wendroff 
equation,  and  a  semi-implicit  4-point  scheme  are  compared. 
The  effect  of  various  boundary  conditions  at  the  right 
boundary  (where  the  differential  equation  requires  no 
boundary  condition  if  n  >  0)  on  stability  is  studied. 
The  non-linear  instability,  when  present,  was  isolated, 
in  each  case,  by  comparison  with  a  calculation  in  which 

-S       2  -s 

4-  ^  was  replaced  by  i^q  "^  ^  '   where  u^   is  the  exact 
solution  of  the  non-linear  equation  (for  the  problems 
chosen,  this  was  known  analytically) .   As  was  expected, 
the  leap-frog  scheme  showed  non=llnear  instability  under 
circumstances  In  which  the  linearized  equation  was  either 
stable  or  weakly  unstable,  but  it  can  be  stabilized  by 
the  introduction  of  a  dissipative  term,  which,  being  of 
the  fourth  order,  does  not  decrease  the  order  of 
accuracy.   The  Lax-Wendroff  equation  was  stable  and  highly 
accurate  in  all  cases.   The  4-point  scheme  was  surprisingly 
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stable,  and  a  tentative  explanation  for  this  Is  given. 

In  part  II,  the  difference  equations  for  coupled 
sound  and  heat  flow.  In  which  the  latter  Is  treated 
implicitly,  are  considered.   The  stability  condition  is 
known  to  be   c  At/Ax  <  1,   where   c   Is  the  Isothermal 
sound  speed.   This  seems  somewhat  paradoxical,  for  if 
the  thermal  conductivity  is  small,  signals  travel  at 
practically  the  adlabatlc  sound  speed  Yj   c,   which  is 
larger  than  c.   For  such  problems,  stability,  as  usually 
defined,  is  inadequate  for  an  actual  calculation,  since 
it  considers  only  the  limit  At,  Ax  ->  Oe   A  practical 
stability  criterion  for  such  problems  is  proposed  and 
Is  applied  to  the  problem  at  hand,  with  complete  success, 
as  judged  by  n-umerical  tests,  which  are  also  described; 
the  criterion  leads,  in  this  problem.,  to  a  stability 
condition  interm.edlate  between  c  At/Ax  <  1   and 
Yj    c  At/ Ax   <   1,   depending  on  various  parameters,  includ- 
ing At  (or  Ax)   itself o   The  condition  is  shown  to  be 
necessary  (by  the  Fourier  analysis)  and  sufficient 
(by  the  energy  method) . 

¥e  wish  to  thank  A.  Wolfe  for  help  with  the  running 
of  these  calculations  and  for  preparing  a  program  for 
finding  the  eigenvalues  of  the  am.plificatlon  matrix  for 
part  II. 
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STABILITY  STUDIES  FOR  DIFFERENCE  EQUATIONS: 

(I)  non-linear  Instability, 

(II)  coupled  sound  and  heat  flowo 

by 
R.  D.  Rlchtmyer  and  K.  W.  Morton 

Part  I:   Non-Linear  Instabilities 

1.   An  example.   Phillips  has  given  an  example  [1]  of 
a  non-linear  partial  differential  equation,  a  difference 
scheme  for  it  that  is  stable  for  the  corresponding  linear- 
ized equation,  and  an  exact  solution  of  the  non-linear 
difference  scheme  which  explodes  as   exp [const ./At|  when 
At  ->  0   for  fixed   t.   A  similar  example  was  given  in 
[2]  and  is  reproduced  here.   The  simple  equation 

(1.1)  ^^+^|-  =   °^        ^  =  u(x,t)  , 

is  considered,  together  with  the  leap-frog  difference 
equation 

/•  n  r^\  n+1    n-1     A  r  /  n   ^  2    /  n   >,  2  -, 

(1.2)  u.    =  u  .    =2  ^^Vl^  "    ^""j-l^    ^       ' 

where  A  =  At/Ax,   and  u.   denotes  u(jAx,nAt),  etc. 

J 

An  exact  solution  of  this  equation  will  be  found,  having 
the  form 
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(1.3)  u^  =   C^  COS  I  j  +  S^  sin  I  j  +  U^  cos  tt  j  +  V5 

it  contains  wave  numbers  7r/(2Z^),  tt/Ax,  and  0.   Substitu- 
tion leads  to  the  recurrence  relations: 

^n+1  ^  ^n-1  ^  2AS^(U^-V)  , 
gn+1  ^  gn-l  ^  ^AC'^CU^-V)  , 
yn+1  __  yn-1  ^  Q 

n 
According  to  the  third  of  these,   U   takes  on  two  constant 

values,  say  A  and  B,   on  alternate  cycles;  then,  the 

first  two  relations  give   the  equation 

(1.4)  C"""^^  -  2  C^  +  C^"^  =   4  A^  (A+V)  (B-V)  C'^ 

n 
and  a  similar  equation  for  S  .   If  the  coefficient  of 

C^  on  the  right  of  (1.4)  lies  between  -4  and  0,   all 

solutions  of  (lo4)  are  bounded;  otherwise  there  is  always 

n  i  n  I 

a  solution  C   such  that   |C  |   grows  exponentially  with 

n,   that  is,  with   t/At. 

To  see  when  the  exponential  solution  can  occur,  note 

that  the  usual  stability  condition  for  the  linearized 

version  of  (1«2)  is 

H  max  [|u|]    <  1   , 

and  we  assume  that  this  condition  is  satisfied;  it  is 
therefore  necessary  that 

A ( I V I  +  max  [|aI,|b|})   <  1   , 
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n 
which  shows  that  the  coefficient  of  C   on  the  right  of 

(1„4)  can  never  be   <  -4 .   If,  furthermore, 
(1.5)  |a|   <   |V|    and    |b|   <   |v|  , 

the  coefficient  can  never  be  >  Oo   Therefore,  if  the 
constant  term  V  of  (1. 3)  is  regarded  as  representing 
a  smooth  solution  on  which  the  other  terms  have  been 
superposed  as  a  perturbation,  then,  if  the  amplitude  of 
the  perturbation  is  below  a  limiting  value  or  threshold 
set  by  (1=5),  it  will  not  grow,  but  if  the  perturbation 
is  initially  above  this  threshold  it  grows  exponentially, 
Whether  this  behavior  is  typical  for  other  forms  of 
perturbation  than  (1.3)  is  unknown. 

2.   The  first  test.   The  mixed  problem 

(2«i)    ^-  +|e^  =  0      O^x^l  ,  t  >  0 

(2.2)      u(x,0)   =   X  0  ^  X  ^  1 

(203)  u(D,t)   =0  t  >  0 

(no  right-hand  boundary  condition  is  needed  for  this 
problem)  is  considered;  its  exact  solution  is 

(204)  u(x,t)   =  x/(H-t)  . 

The  corresponding  linear  problem,  in  which  (2.1)  is 
replaced  by 
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(2.5) 


St   ^  +  T+t^  ^  =  0 


has  the  same  solution.   The  leap-frog  and  Lax-Wendroff 
conservation  systems  for  (2.1)  are,  respectively 
(?\  =  At/Ax,   as  in  Section  1), 


(2.6) 


n+1 

u  . 
J 


u 


n-1   A 


A  r  /  n   ^2    ,  n   ,2-1 


(2.7) 


n   ?\  r  ,  +,  2 
u.  -  2  [  (w  ) 

where 

1  /  n    ,   n, 

2  ^Vl  ""^-^ 

and 

1  /  n  ,   n   , 

2  ^^-  -^  ^--1^ 


(w")2]   , 


A  r  /  n   ,2    ,  nv  2-, 
^  t(u^.^^)   -  (u.)  ] 


A  r  /,  n,  2    ,  n   ,2. 
2  [  (u.)   -  (u.  ^)  ] 


The  leap-frog  and  Lax-Wendroff  linearized  schemes  for  (2.5) 
are,  respectively. 


(2.8) 


n+1 

u  . 
J 


u 


n-1 
J. 


n 


Ax   /  XI 


""-l'  ' 


(2.9) 


n 
u  .  - 
J 


7\  X 


+ 


Hht+  -|  At 


(w  -  w  ) 


where 


X  +  -o  Ax 
2  ^^j+1  +  ^j)  "  >^  —T+t ^^-+1  -  ^j^ 


1  /  n 


and 


X  -  -?r  Ax 
2  ^^-  +  Vl^  "  ^  -T+t ^^-  -  ^^J-l^  ' 


1  /-  n  .   n 


here  x  =  jAx  and  t  =  nAt. 

In  the  runs  reported  below,  A  =  1,   so  that  the 
stability  condition 


At 

-T—  max 

Ax 


{i"il 


is  satisfied  for  all   t  >  0. 

For  these  difference  equations,  a  boundary  condition 
is  unfortunately  needed  at  x  =  1,   i.e.  at   J  =  J,   where 
Ax  =  l/j.   In  the  first  test,  we  simply  set  u,   equal 
to  the  value  of  the  exact  solution  (2.4)  at  x  =  1, 
t  =  (n+l)At. 

To  study  the  growth  of  perturbations,  the  initial 
values  were  computed  by  the  formula 

u°  =  jAx(l+0(e  .  -  h)    , 


(2,10) 


u"-'-  =   u^  /  (1-At)  , 


where  ^^,...,1-^  are  random  numbers  uniformly  distributed 
J-    '  J 

between  0  and  1,  and  9  is  a  sm.all  parameter  that  can  be 
varied  from  one  run  to  another. 

In  Figure  1  is  plotted  the  average  of  the  perturbation, 
that  is, 

^^  f '^j  "   ift'    0  <  j  <  J,  X  =  jAx,  t  =  nAtj 

as  a  f\j.nction  of  the  cycle  n^umber  n,   for  the  four 
difference  schemes,  for  9  =  0.3,  (the  mean  initial 
perturbation  therefore  was  7  -p  /o)  .   It  is  seen  that 
(a)   the  non-linear  leap-frog  equation  explodes  abruptly 
after  cycle  200,   (b)   in  the  linearized  leap-frog 
equation,  the  perturbation  neither  grows  nor  decays, 
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(c)   the  perturbation  is  damped  out  by  the  Lax-Wendroff 
equation,  especially  in  the  conservation  form;   (b)  and 
(c)  are  in  accord  with  the  fact  that  the  amplification 
factors   g   (see  [3],  p.  55)  lie  on  the  unit  circle  in 
the  complex  g  plane  for  the  leap-frog  scheme  and  inside 
it  (except  at  the  point  g  =  1)   for  the  Lax-Wendroff 
equations,  which  are  therefore  dlssipative,  in  the  sense 
of  Kreiss  [5] . 

The  abruptness  of  the  explosion  for  the  top  curve 
(a  pure  exponential  rise  would  be  represented  by  a  :v:;.>;. 
straight  line  in  the  graph)  suggests  that  the  threshold 
phenomenon  noticed  in  the  example  given  in  Section  1  may 
be  typical,  for  during  the  first  I50  cycles  the  average 
perturbation  did  not  change,  while  the  magnitude  of  the 
exact  solution  x/(l+t)   decreased  by  a  factor  6; 
presumably  the  relative  amplitude  of  the  perturbation 
(i.e.  relative  to  the  smooth  part)  exceeded  the  threshold 
around  cycle  l80-200. 

Runs  with  other  values  of  the  parameters  0,r  were 
similar. 

3-   Influence  of  boundary  conditions  on  the  non-linear 
leap-frog  equations.   At  the  suggestion  of  Kreiss,  who 
pointed  out  that  the  right  boundary  condition 
(Uj  =  x/(l+t))   could  be  Itself  a  cause  of  instability, 
a  second  test  calculation  (modification  3  in  the  series) 
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was  made,  in  which  the  non-linear  leap-frog  equation  was 
solved  in  the  same  manner  as  in  the  first  test  except  that 
5  different  choices  of  right  boundary  condition  were  used, 
as  follows: 

1.  u^''"''"  =  ulo  (horizontal  extrapolation,  order  0), 

2.  Uj"*"-^  =   2  Uj^2  "  ^jlj         (horiz.  extr . ,  order  1), 


3.  u^+1  =   2u%-u^-l  (115°  extrap.),   /^^ 

"4.  uf^     =  I  t^^r^)^  -  ^^5:2)^]  ("explicit"), 

n+X 
5.  U-,    is  taken  as  the  solution  of  the  equation 

n+1    n-1  ,  A  r,  n+lv2  ,  n+lv2-,      „  mi-,„„t  ,- „^. +-im 
u^   ~  ^j   +"2  '- '  J  '    "^  J-2^  J   =   0  (  implicit  j. 

The  last  two  are  difference  approximations  to  the 
differential  equation.   The  apparently  strange  choices 
of  the  indices  were  made  because,  in  the  leap-frog  scheme, 
the  set  of  values  u .   with  n+j   even  are  completely 
uncoupled  to  the  set  with  n+J   odd;   that  is,  there  are 
tivo  completely  independent  calculations  in  progress 
simultaneously,  and  it  was  felt  that  to  couple  them 
together  by  the  boundary  condition  could  only  introduce 
unwanted  complication  ;    hence  the  above  equations  were 
chosen  to  preserve  parity  of  superscript  plus  subscript 
in  each  equation.   (One  of  the  two  independent  calcula- 
tions could  have  been  suppressed,  of  course,  at  the  cost 
of  a  slight  bit  of  additional  programming.) 
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Figure  2  shows  that  all  these  calculations,  except 
number  5,  with  the  ^5   extrapolation  for  boundary  condition, 
explode,  some  sooner  than  the  others .   In  a  second  run, 
with  a  larger  initial  perturbation,  method  niomber  3  also 
exploded. 

It  is  possible,  of  course,  that  some  other  boundary 
condition  than  any  of  those  tried  might  have  a  stabiliz- 
ing effect. 

4 .   Boundary  conditions:   non-linear  Lax-Wendroff  equation. 
In  view  of  the  violence  of  the  explosions  reported  above, 
a  third  test  was  made  (mod.  4)  in  which  the  Lax-Wendroff 
method  (2.7)  was  subjected  to  the  same  five  boundary 
conditions  described  above.   All  were  stable.   With 
A  =  1  and  0  =  0.5   (a  fairly  large  initial  perturbation), 
the  perturbation  had  dropped  at  cycle  300  to  4*10  -  for 

the  first  method  (zero-order  horizontal  extrapolation)  and 

r        -4 
to  6 '10   for  the  others. 


5'   Boundary  conditions:   linearized  leap-frog.   In  this 
test  (mod.  5),  the  linearized  leap-frog  equation  (2,8)  was 
subjected  to  the  same  five  boundary  conditions  (Section  3) 
at  X  =  1;  the  average  perturbation  is  plotted  against  n 
in  Figure  3.   In  no  case  was  there  an  exponential  growth, 
but  in  the  first  two  methods  there  was  a  slow  growth  of 
the  perturbation,  which  rendered  the  computed  values 
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meaningless  after  about  50  cycles,  and  it  is  possible 
that  this  would  have  happened  in  the  other  three  methods 
also.  If  they  had  been  continued  longer.   Evidently,  the 
best  that  one  can  hope  for  In  the  leap-frog  method  Is 
that  errors  will  neither  grow  nor  decay,  and  this  Is 
presumably  because  the  amplification  factors  all  lie  on 
the  unit  circle  in  the  complex  plane. 

The  growth  Is  shown  in  a  log-log  plot  in  Figure  4, 
for  boundary  conditions  1  and  2;    the  graphs  suggest 
that  the  growth  is  asymptotically  with  a  low  power 
(~/l.O  and  0.2  respectively)  of  n. 

6.   A  four-point  scheme.   In  order  to  be  completely  rid 
of  the  effects  of  a  boundary  condition  at  x  =  1,   the 
Initial- value  problem  (2.1),  (2.2),  (2.3)  was  next  solved, 
in  "modiflcatlDn  6, "  by  the  scheme 

,r   ,,    n+1  ,   n+1    n     ,  n  .  A  r,  n+1^2  ,  /  n   ,2 
(6.1)   u^._^-j_  +  u.    -  u^.^^  -  ^j  +^  f  (^j+l)   +  (^j+l) 

,  n+1, 2    /  n,2i      „ 
-  (Uj   )   -  (Uj)  ]   =  0   . 

n+1 
If,  at  each  cycle,  the  unknowns  u ,    are  solved  for 

J 

in  order  of  Increasing  j,  starting  with   j  =  1 

(ui   =0  by  (2.5)),   the  equation  is  effectively  explicit, 

n+1 
in  that  it  then  contains  only  one  unknown,  namely  u.^  , 

for  each   j,   and  one  merely  has  to  solve  a  quadratic 

n+1 
equation  to  find  '^■^ii    >       (3   =  l,.o.,J). 
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Equation  (6.1)  Is  similar  to  the  leap-frog  equation 
In  that  the  amplification  factors  all  lie  on  the  unit 
circle,  for  the  corresponding  linear  equation  with 
constant  coefficients. 

A  dlsslpatlve  term,  with  an  adjustable  coefficient, 
was  also  added  to  (6.1),  and  It  was  our  Intention  to 
Investigate  the  smallest  value  of  this  coefficient  that 
would  lead  to  stability.  Much  to  our  surprise,  however, 
the  procedure  turned  out  to  be  completely  stable  without 
any  dissipation  at  all.   For  0  =  0,3,   A  =  1.0,   J  =  30, 
the  mean  absolute  perturbation  decreased  from  an  initial 
value  of  7.5  •  10"^  to  between  2.6-10"^  and  2.8'>10"  , 
at  n  =  300,   depending  on  the  coefficient  of  the  dlsslpa- 
tive  term;  the  lowest  value,  2.61 '10    was  for  the  case 
of  no  dissipation. 

A  tentative  explanation  for  the  stability  of  this 
scheme  (later  rejected  however  —  see  Section  7  below) 
was  arrived  at  by  a  closer  examination  of  the  amplifica- 
tion factor,  as  follows:   to  linearize  (6.1),  we  write 
u.  +  e.   in  place  of  u.   in  (6.1),  where  now  u. 
stands  for  a  smooth  solution  of  (6.1)  and  b^     for  a 
perturbation.   The  contribution  of  the  square  bracket 
term  to  the  first  variation  of  (6.1)  is  the  quantity 

r c   ns  V     T,  ^,,^+1  ^+1   ,  ,  ri    n       n+1  n+1   n  n^ 

(6.2)     X  =   Muj^-L  ej^^   +Vl'j+1   ~^J   'j   "VJ^ 

We  approximate  u.^   as 
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n+  - 
n+1  r^      I       ,  At  Su  ,  Ax  Su,    2 

and  similarly  for  the  other  terms  In  (6.2);  then, 

,.  a.  ^  /  n+1  ,   n      n+1    r\x 
X     ^     Au(B  .^,  +  s  .^,  -  s  .    -  e  .) 

,  A  .,  Su  ,  n+1    n      n+1  ,   nN 
+  2  ^t^  ^^-+1  -  ^J+1  -  ^j    -*-^J^ 

+  ^  Ax#i  (s^t!  +  s^^,  +  B^-^1  +  e^)  . 
2     dx    j+1     j+1     J       J 

This  X,  together  with  contributions  from  the  first  four 
terms  of  (6.1),  must  be  equated  to  zero.   To  find  the 
amplification  factor  g  =  g(k)   for  the  Fourier  component 
e^^,   we  substitute   (g)'^(  e^^^^)  •^"   into  the  equation  of 
variation;  this  gives 

g[  (1  +  A  ^  1^)  cos  a  +  iA(u  +  ^  |^)sin  a] 


(6.3) 


=   (1  „  A  -f  1^)  cos  a  =  iMu  -  ^  |H)sin  a  =  0, 


where  a  =  kAx/2.   This  equation  is  of  the  form  gA  =  B, 
and  we  will  have   |g|  <  1   if   |b|  <  |a|,   or  if 
|a|^-  |b|^>0.   From  (6-3) , 

I  n  I  2    I  -n  I  2     o•^  A   Su  ,  o^  2   . ,  ^u 
A    -   B     =   2A  Ax  -r—  +  2A  u  At  -jr-  , 
'  '     '  '  ox  dt 

or 

(6.4)      lAl^  .  |B|2   =   2  At  (|H  +  A^  u|H)  , 

For  the  case  at  hand,   u  =  x/(l+t),   and  {6  A)    is  positive; 

i.e.,  the  values  of  g  all  lie  inside  the  unit  circle 

except  for  the  point  g  =  1. 
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7.  The  four-point  scheme  with  dissipation  and  a  forcing 
term.  To  find  a  solution  in  which  {6  A)  is  negative,  at 
least  for  some  x  and  some  t  (our  objective  is  still 
to  find  a  scheme  that  is  marginally  stable  in  the  linear- 
ized form  but  subject  to  non-linear  instabilities,  which 
can  be  made  stable  by  adding  a  suitable  dissipative  term 
--  see  next  section  for  a  successful  attainment  of  this 
objective),  we  considered  the  equation 

(^•1)  It-  ^^4  =  ^ ' 

where  f   =  ■i;i/(x,t)   is  a  so-called  forcing  term. 

If  ?//  =  -^  sin  2irx,      then  (7.1)  has  a  stationary 
solution  u(x,t)  =  sin  irx.      In  modification  7  we  took 
this  expression  for  Tp     ^-i^d  we  took  u(x,0)  =  sin  irx     as 
initial  condition.   Then  -^p  =  0   and,  by  {6oh),       \k\   -|b| 
is  positive  for  0  <  x  <  1/2  and  negative  for   1/2  <  x  <  1. 
The  von  Neumann  condition   |g|  _<  1  +  0(At)   is  satisfied 
everywhere,  but,  according  to  the  foregoing  considerations, 
one  might  expect  non-linear  instabilities  to  develop  in 
the  region  1/2  <  x  <  1,   where   |g|   actually  exceeds  1 
by  0(At). 

The  difference  equation  was  identical  with  (6.1), 

v  1 

except  that   2At(-p  sin  27r(j  +  -p)Ax)   was  put  on  the 

right  side,  instead  of  0.   To  our  astonishment,  this 

procedure  was  also  stable  without  any  dissipation; 

the  average  perturbation  after  500  cycles  was  ^   4  "10"^ 

with  no  dissipation  and  up  to   4.5-10"^  with  the 

maximum  possible  amount  of  dissipation.    These 
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figures  are  about  I5  times  larger  than  those  of  the 
preceding  section,  possibly  indicating  a  slight 
destablizing  effect  of  the  region  1/2  <  x  <  1,   where 

|g|  >   1- 

It  is  conjectured  that  the  stability  of  this 
scheme  can  be  explained  as  follows:   the  equation  can 
be  written  as 

/  „  ^  X  Su  ,    Su      , 

where   c  =  c(u),   the  velocity  of  propagation  of 
perturbations,  is  equal  to  u   itself,  and  hence  is 
positive.   Therefore,  perturbations  arising  in  the 
region   1/2  <  x  <  1   are  swept  to  the  right  boundary, 
where  they  disappear,  before  they  can  be  amplified 
by  more  than  a  small  factor,  so  that  the  threshold  for 
the  non-linear  effects  is  never  exceeded.   This  conjec- 
ture is  somewhat  supported  by  the  results  of  the  next 
test  (Section  8,  below),  where  strictly  periodic  boundary 
conditions  are  used,  so  that  perturbations  are  kept  in 
the  system  forever;  then,  some  positive  dissipation  is 
indeed  required  to  prevent  non-linear  explosive  instabil- 
ity. 

The  tendency  for  perturbations  to  be  swept  to  the 
right  boundary  and  disappear  was  also  observed  by  Parter 
[6].   For  the  Lax-Wendroff  scheme  applied  to  an  initial- 
value  problem  based  on  (7.2),  he  found  that  the  equations 
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behaved  in  a  stable  manner  and  gave  accurate  results 
even  with   c  At/Ax  slightly  greater  than  1   (I.O5,  in 
his  example),  provided  the  number  J   of  net-points 
was  not  too  large  (not  over  about  5O)  .   Parter  called 
this  pehnomenon  "pseudo-stability." 

8.   Leap-frog  with  dissipation,  forcing,  and  periodicity. 
Next,  in  "modification  8,"  the  equation 

(8.0)  Jt  ^  +^^  =     TT  sin  47rx 
was  approximated  as 

(8.1)  u^+^  =  u^-^  -^  [(u^   )2-(u^   )2]  +  (tt  sin  47rJAx)2At 

e  ,  n-1    ),  n-1  ,  ^   n-1    ),  n-1  ,   n-lx 
-  Tf  (^3+2  -  "^^j+l  +  ^^j    -  ^^j-1  +  ^J-2)  • 

In  place  of  boundary  conditions,  a  periodicity  condition 

was  used:   in  each  cycle,  (8.1)  was  used  for   j  =  2,3,..., 

T ,  n    ^i„     n+1     ,   n+1         ,      ^  ,     n+1     -, 
J+1;   then  u^  and  u^     were  set  equal  to  u^    and 

n+1   „  -,   n+1     -,   n+1         ,      ,  ,     n+1     -, 
u^  -,,   and  u,  p  and  u^  ^  were  set  equal  to  Up    and 

n+1 
u^   . 

The  last  term  in  (8.1),  which  approximates 

=  (s/4)  (S  u/Sx  )(Z\x)  ,   has  a  dissipative  effect;  it  makes 

the  amplification  factors  lie  inside  the  unit  circle,  if 

A  and  e   are  property  chosen.   It  will  appear  that  the 

dissipative  effect  of  this  term  is  achieved  at  the 

expense  of  a  slight  decrease  of  the  allowable  value  of  A, 

just  as,  for  example,  when  an  artificial  viscosity  term. 
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is  used  for  the  calculation  of  shocks-   (See  [5]j.  P-  2l8 
ff .) 

The  stability  analysis  will  be  given  first  for  the 
case  of  no  dissipation,   e  =  0.   Equation  (8.1)  is  first 
linearized  by  writing  u  .  +  t]  .  in -place  of  u.,   where 

J         J  J 

u.   is  a  smooth  solution  of  (8.1)  and  ri  .   is  a  perturba- 

n 
tion.   To  first  order,  the  equation  for  r\  .      is 

n+1    n-1  ,  , ,  n   n      n   n   ,      „ 
ri.    -   r\  .         +  7\[-u.  .  ,  .T]  .  ,  ^    -   u  .    -.r[  .    ^  I      =   0. 

Since  u'? , -,  f^  (u  +  Ax -r^)  .,   the  above  is  approximately 
j+1  ox  J 


n+1      n-l,-,,n  n      ^     ,  ^    .^6u    ,n      ,n      -,  ^ 

■n  .      -ri  .         +  Au  ( ri  .  ,  T  -Ti  .    -,  )    +  A  Ax  -r-    HI  .  ,  i  +11  .    n         =      0  . 
',1  '.1  ',1+1     ',1-1'  dx       'j+1     'j'-l 


To  find  the  amplification  factor  g,   substitute 
(g)  (e    )'^      into  this  last  equation;  this  gives 

(8.2)     g^-1  +  (2iAu  sin  ct  +  2A  Ax -^  cos  a)g  =  0  , 

where  a  =  kAx.   For   I  Au  I  <  1  and  -r—  =  0,   the  roots 

'  —         ox 

of  (8.2)  lie  on  the  unit  circle  on  either  side  of  the 
imaginary  axis  and  their  sum  is  pure  imaginary  (=  -2iAu- 
sin  a) .   The  effect  of  the  extra  term  in  Su/Sx   is  to 
shift  the  roots  slightly  laterally  so  that  one  lies 
slightly  outside  and  the  other  slightly  inside  the 
unit  circle,  regardless  of  the  sign  of  bu./bx.      Since  the 
shift  is  proportional  to  A  Ax  =  At,   the  von  Neumann 
condition   | g |  _<  1  +  0(At)   is  satisfied,  but  non-linear 
instability  may  nevertheless  be  expected. 
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The  effect  of  the  dlsslpatlve  term  in  (8.1)  will 
now  be  considered;  for  simplicity  we  set   Su/Sx  =  0. 
Then,  (8.2)  is  replaced  by 

P  2 

(8.5)     g   -  2  i  Au(sin  a)g  +  (l-e(l-cos  a)  )  =   0  . 

Since  the  product  of  the  roots  must  be  <  1  for  stability, 
we  assume  that  0   <_  e    <■   1/2   in  any  case.   When  the 
discriminant  D  of  this  equation,  given  by 

(8.4)     ^  =   -  A^u^  sin^  a  +  1  -  e  (1-cos  a)^  , 

is  positive,  the  roots  of  (8.3)  are  symmetrically  located 
relative  to  the  imaginary  axis,  and  each  has  magnitude 
lA-£(  1-cos  a)   .    As   a  varies  from  0  to  2Tr,   g  traces 
out  a  path  in  the  complex  plane  consisting  of  one  or 
more  symmetric  curves  lying  inside  the  unit  circle 
together  with  one  or  more  line  segments  lying  on  the 
imaginary  axis.   These  segments  may  stick  out  of  the 
unit  circle;  for  example,  if  Au  =  1,   a  =  ir/2,   e  =  1/2, 
we  find  g  =  1(1  +-yi/2)  ;   to  prevent  this,   Au  must 
be  kept  somewhat  below  the  value  1  that  applies  when 
e  =  0. 

Since  our  objective  is  to  provide  dissipation,  it 
is  not  enough  merely  to  prevent  these  segments  from, 
extending  outside  the  unit  circle  --  we  want  the  values 
of  g   to  lie  well  inside  it,  except  for  a  in  the 
neighborhood  of  0  or  ir .      Consideration  will  therefore 
be  restricted  to  values  of  e   and   Au  such  that  the 
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discriminant  D  of  (8.3)  is  positive,  so  that  the  path 
traced  out  by  g   consists  only  of  the  symmetric  curves. 
We  write  D/4  =  1  -  f(cos  a),   so  that 

(8.5)      f(p.)   =   £(l-M-)^+A^u2  (1-H^)  , 

and  we  require  that   f(|a)  _<  1   for   -1  _f  M-  _^  !•   This 
leads  (we  omit  details)  to  the  restrictions 

0  _<   e  ^  1/4 

(8.6) 

A^u^  <     (1  +yT^^)/2  . 

The  last  test  in  this  series,  mod.  8,  was  made  to 
find  out  for  what  values  of  e   and   A   (max  {u}  =  1 
for  the  exact  solution),  within  the  restriction  (8.6), 
the  calculation  would  be  stable  against  the  non-linear 
effects.   The  initial  function  was  taken  to  be  u(x,0)= 
sin  27rx,   plus  a  small  perturbation,  and  the  solution 
of  (8.0)  is  then  the  stationary  solution  u(x,t)  =  sin  27rx. 
Equation  (8.1)  was  solved  for  a  few  hundred  cycles,  for 
various  values  of  e   and   t. 

For  example,  with  A  =  O.9,   (8.6)  holds  if 
0  _<  £  ^  approximately  0.125;  in  a  run  of  3OO  cycles  the 
calculation  was  stable  for  e   =  .02, .05, .04, .06, .08, .10, 
.125, .15;   the  average  perturbation  dropped  from  To 7° 
to  2  or  3^0.   For  e  =  0  and  .01,  the  calculation  exploded, 
that  is,  values  of   |u|  >  1000  were  obtained  and  the 
calculation  automatically  terminated.   With  no  dissipation 
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(e  =  0),   the  explosion  took  place  between  cycles 
n  =  50   and  n  =  100;   with  e   =  .01,   It  took  place 
between  n  =  I5O  and  n  =  200.   Analogous  results 
were  obtained  for   A  =  0.8  ,   for  which  the  allowed 
range  of  e   by  (8.6)  is   0  _<  e  _<  0.2,   and  for 
A  =  0.5,   for  which  it  is  the  full  range  0  _<  e  j<  1/4. 
It  appears  that,  at  least  for  the  simple  equation 
(8.0),  the  leap-frog  scheme  can  be  stabilized  by  a 
small  perturbation.   (It  is  a  small  quantity  of  fourth 
order,  whereas  the  error  of  the  first  two  terms  of 
(8.1)  is  of  third  order.) 


Part  II.   Coupled  Sound  on  Heat  Plow. 

9.   Statement  of  the  problem.   The  equations  for  small- 
amplitude  or  acoustic  motion  of  a  gas  whose  thermal 
conductivity  is  c      are 

o 
Se     ^  S  e     Su 

^  =  ^—2  -  c  ^  ; 

dx 
here,   7   is  a  constant  >  1   and   c   is  the  isothermal 
sound  speed.   For   o-  =  0,   the  system  (9«1)  reduces  to 
the  wave  equation  with  propagation  speed  equal  to  Yyc, 
the  so-called  adiabatic  sound  speed.   (See  [5],  p.  170  ff. 
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--  but  note  the  error  in  sign  in  the  first  equation  (8.8) 
and  the  first  equation  (8.9).) 

In  [^] ,    it  was  conjectured  that  the  difference 
approximation  to  (9-1)  given  by  the  equations 

n+1   n        n       n      /   -,  %  /  n       n     ^ 
^J   -^.1   _   ,  ^-^  1/2-^1-  1/2-^^-^^  ^^1+  1/2-^ j-  1/2^ 
At    "   °  Ax 

n+1     n  n+1    n+1 

(q    p]  J  +1/2   ,1+  1/2   _   p   J+1     J 

^^'^'  At  ~   ^      Ax 

n+1     n  n+1   _„  n+1     n+1 

^j+  l/2'^j+  1/2  _     ^J+  V2~^^J+  1/2+^j-  1/2 

At -   "^ —^ 

Ax 

n+1   n+1 
u  .  ,  -,  -u  . 

_  p   J+1   J 

^  Ax 

of  which  the  first  two  are  effectively  explicit  while 
the  third  is  implicit,  has  as  its  stability  condition 
c  At/Ax  <  1.   The  correctness  of  this  conjecture  was 
confirmed  in  I962  by  Morimoto  [4]  .   In  the  limit  <T ->   0, 
however,  the  stability  condition  is  the  more  stringent 
inequality  yyc  At/Ax  <  1,   and  this  suggests  that, 
unless  Ax  <<  (T/c,      there  should  be  a  practical 
stability  condition,  for  a  finite  net,  involving  the 
two  dimensionless  constants 

c  At  o-At 

^^  (Ax) 

and  7;      it  would  presumably  be  intermediate  between 

the  two  conditions  above,  namely,  it  would  be  of  the 

form  V  <  Vq,   where   l///  <   Vq   <   '^,      where   Vq   is 
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some  function  of  7   and  [i . 

10.   Modified  von  Neumann  condition.   In  [3],  It  was 
pointed  out  (p.  1^4  and  p.  I88)  that,  for  some  problems, 
stability  Is  not  enough  for  practical  utility  of  a 
method.   Although  stability  implies  convergence  in  the 
limit  At  ->  0,   it  can  happen  that,  for  the  finite  At 
used  in  practice,  errors  are  nevertheless  unacceptably 
amplified.   If,  in  the  sound  and  heat  flow  problem,   (T* 
is  so  small  that  it  is  impractical  to  choose  Ax  <<  cr/c, 
or  even  Ax  ~  (T^/c,      and  we  are  forced  to  accept  Ax  >>cr/c, 
then  the  system  (9'2)  must  behave  almost  exactly  like 
the  wave  equation  with  propagation  speed  ■Yj   c;   there- 
fore, if   cAt/Ax   exceeds   I//7  t)y  an  appreciable  amount, 
one  must  expect  errors  to  amplify  just  as  they  do  when 
the  Courant-Prledrlchs-Lewy  condition  is  violated  for 
the  wave  equation. 

The  remedy  indicated  in  [5]  Is  to  require  that  no 
component  of  the  approximation  be  allowed  to  grow  more 
rapidly  than  the  most  rapid  possible  growth  of  an  exact 
solution.   We  therefore  require  that  the  constant  implied 
by  the  term  0(At)   in  the  von  Neumann  condition 
|g|  jf  1  +  0(At)   be  somehow  restricted  by  physical 
considerations . 

A  possible  formalization  of  this  idea  is  the  follow- 
ing.  We  consider  a  pure  initial-value  problem  of  p 
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partial  differential  equations  of  first  order  in  t 

->  ->> 

for  p   functions  of   t  and   d   space  variables,   u(x,t), 

■*■  -> 

where  u  has  p   components  and  x  has   d  components. 

According  to  von  Neumann's  procedure,  the  equations  are 

linearized  and  the  coefficients  are  treated  as  locally 

constant;  then 

It  ^  (^'^)   =  ^^^  '■'"   ^^  ^  ^^'^)  ' 

where  D(q)   or  D(q-,  ,  .  .  . ,  q  , )   is   apxp  matrix  whose 
elements  are  polynomials  in  q  ,  ...,q, .   After  the  Fourier 
trans  format  ion 

u  (x,t)   =   /  d  k  v(k,t)  e    '     , 

the  formal  solution  of  the  initial-value  problem  (see 
[3],    p.  52)  is 

V  (k,t)   =   e      '  v  (k,0)  ; 

therefore,  if  we  call 

C  \  ^  ^2 

(10.1)    m  =  max  j Re  A   A  =  e-val.  of  D(ik),  all  real  kj, 

then  the  solution  cannot  grow  faster  than   e   ,   possibly 
multiplied  by  a  polynomial  in  t.   We  therefore  take  as 
a  modified  von  Neumann  condition 

|g|   =   |g(At,k) I  ^   e    , 

where  g   is  any  eigenvalue  of  the  amplification  matrix 
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of  the  difference  scheme  under  consideration,  and  m 
is  given  by  (10 .1) . 

In  [3],  an  approximate  application  of  this  condition 
was  made  (p.  I35)  to  a  problem  of  neutronics,  in  which 
m  >  0,   and  a  strict  application  was  made  to  a  problem 
of  a  vibrating  bar  under  tension,  where  m  =  0   (p.  I88) . 


11,   Amplification  factors  for  the  sound  and  heat  flow 
problem.   To  find  the  eigenvalues   g   of  the  amplifica- 
tion matrix  (which  we  shall  not  actually  write  out),  we 
substitute 


u 


0 


,  ,n  ,  ikAx^J 
(g)   (e    )^ 


w 


0 


'0 


where  u^,    w„  and  e^     are  constants,  into  (9.2).   For 
a  non-trivial  solution  Up^,  w„,  e„   to  exist,  we  must 


have 


(11.1) 


g-1 


-2iv  sin  a   2iv(7-l)sin  a 


-2igvsin  a  g~l 


2igvsin  a,     0 


0 


g-l+2[i(l-cos  2a)g 


=  0 


Since  sound  waves  do  not  grow  in  amplitude  (in  fact, 
they  are  damped  out  by  the  action  of  the  heat  flow),  the 
quantity  m  defined  by  (10.1)  is  zero  in  this  problem, 
and  we  require  that   v   be  such  that   |  g  |  _<  1   for  all 
solutions  of  (11.1),  for  0  <  a  <  27r. 


-  26  - 


Figure  5  shows  the  locus  In  the  complex  g  plane 
of  the  roots  of  equation  (11.1)  as   a  varies  from  0 
to  2t,      for  particular  values  of  v,    \i,    and  y ;      it 
consists  of  a  curve  and  two  segments  of  the  real  axis. 
This  case  was  very  close  to  the  stability  limit,  and 
the  left-hand  segment  of  the  real  axis  extends  very 
slightly  outside  the  unit  circle.   The  left  endpoint 
of  this  segment  is  one  of  the  roots  of  (11.1)  for 
a  =  Tr/2;   it  moves  to  the  left  as   v   increases.   Assum- 
ing that  this  is  the  case  also  for  other  values  of  v,  pi, 
and  -y ,      we  can  find  the  m.odified  von  Neumann  condition 
by  substituting  -1  for  g  and  1  for  sin  a     into  (11.1) 
to  obtain  a  limiting  value  of  v.   This  will  in  any  case 
give  us  a  necessary  condition  for  stability  in  the 
modified  sense,  and  a  separate  argum.ent,  given  in 
Section  13  below,  shows  that  this  condition  is  also 
sufficient.   The  substitution  gives 

-1  -  2|j.  +  v'^7  +  2  (j.v^  =      0   , 

hence,  we  take  as  the  practical  stability  condition 


(11.2)         V   <  /^^ 


7+2[i 

As  was  expected,  stability  is  determined  by  the 
isotherm.al  sound  speed  for  [x  »  1   (large  o"  ) ,  by  the 
adiabatic  sound  speed  for  [j,  «  1,   and  by  an  inter- 
mediate speed  for  |-L  ~  1 .   This  can  be  correlated  to 
some  extent  with  the  behavior  of  actual  sound  waves  of 
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wavelength  comparable  with   2Ax.   Substitution  of 
u^    exp   |'l(kx-a3t)?   for  u  and  similar  expressions 
for  w  and   e   Into  the  differential  equations  (9'1) 
yields  the  equation 

(11.3)  1^  (^^-1)  +  M<l'^-7)   =   0 


as  the  condition  for  a  non-trlvlal  solution;  here, 

,-,-.    2|V  (f,     -     iL.     -     wave   speed 

^       '    '  ''^  ~   kc  ~   Isothermal  sound  speed 

For  given  real  wave  number   k   (the  wavelength  is   27r/k) , 
the  three  roots  of  (11. 3)  represent  a  sound  wave  moving 
to  the  left,  a  sound  wave  moving  to  the  right,  and  a 
stationary  sinusoidal  temperature  distribution  that  is 
slowly  decaying  by  heat  conduction.   For  k  c/c  »   1, 
that  is,  for  wavelength   <<  c/cr  ,   (|)  %  +1   for  the 
roots  of  (11.3)  that  correspond  to  sound  waves;  the 
waves  therefore  travel  with  speed   c.   For  k  o/c   <<  1, 
(j)  ^  +  Yj    ;      here,  the  waves  travel  with  speed  -/y  c. 
For  k  (T/c  ~  1,   the  waves  have  an  intermediate  speed 
and  are  damped,  because   (j)   then  has  a  rather  large 
negative  imaginary  part.   To  illustrate,  with   k  =  tt/Ax, 
corresponding  to  the  shortest-wavelength  disturbance 
that  can  appear  in  a  calculation,  if  we  also  take 
(T'/(cAx)  =  1,   7  =  3^   solution  of  (11. 3)  gives 

(j)  =  +  1.079  -  .35091,     0  -  2.441  . 
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These  results  are  in  qualitative  accord  with  the 
Idea  that  stability  Is  Influenced  by  the  behavior  of 
short-wavelength  disturbances. 

12.   The  numerical  tests.   The  system  (9.2)  was  solved, 
for  given  7,  v,    \x,    J,    with  initial  and  boundary  data 
given  by 


(  "  \ 


w 


l^  / 


0 
0 
0 


for  J  <  j/3,  n  =  0, 

and  J  =  0  ,  all  n, 

for  j  >  J/5,  n  =  0, 

and  J  =  J  ,  all  n. 


If  C    were  =  0,  the  corresponding  exact  solution 
of  the  differential  equation  would  be  a  step  wave,  or 
infinitesimal  shock  wave,  moving  to  the  right  with 
speed  7/7  c.   With  (T'  /  0,   we  do  not  know  the  exact 
solution,  but  an  approximate  solution  shows  that  the 
wave  still  moves  to  the  right  with  speed  -/y    c,   but 
becomes  gradually  smoothed  out;  it  has  approximately 
the  shape  of  the  normal  distribution  function 
I        exp  5-y  /2w  j  dy   with  a  width   w  =  Vcr  t . 

The  calculation  was  so  programmed  that  whenever 
enough  time  had  elapsed  so  that  the  wave  ought  to  be 
at   J  =  2j/5,   the  wave  was  shifted  bodily  to  the  left 


-  29 


by  one  third  of  the  distance   JAx  --  that  is,  values 
for  0  <  j  _<  2j/3  were  replaced  by  those  for 
J/3  <  J  ^  J 7   and  values  for  2J/3  <  j  _<  J  were  set 
to  zero  —  in  this  way,  the  calculation  could  be 
continued  indefinitely »   Figures  6  and  7  show  typical 
velocity  profiles .   The  speed  with  which  the  inter- 
section of  the  curve  with  a  horizontal  line  u  =  O.5 
moved  to  the  right  agrees  with  ^/y  c   to  about  0.4 /b. 
The  agreement  between  the  two  curves  in  Figure  7  (they 
are  for  the  same   t  but  different  At)   suggests  that 
the  calculation  is  fairly  accurate. 

Figure  8  summarizes  the  results  of  a  set  of  calcula- 
tions for  7  =  3.0.   Each  solid  dot  indicates  a  run  that 
was  stable  in  the  sense  of  giving  smooth  curves  like 
those  of  Figures  6  and  7  after  the  wave  had  run  a 
distance   5/5   times   JAx,   and  each  cross  indicates  a 
run  that  became  unstable  in  the  sense  that  an  oscilla- 
tion of  amplitude  %  2.0  times  the  average  value  of  u 
appeared,  as  in  Figure  S,    at  which  time  the  calculation 
was  stopped;  in  each  case  this  happened  in  the  first 
five  cycles  and  the  oscillation  always  appeared  near 
the  front  of  the  wave. 

The  theoretical  stability  limit  given  by  (11.2) 
is  also  shown  in  Figure  8  and  is  seen  to  be  in  agreem.ent 
with  the  numerical  experiments.   Two  calculations  that 
were  exactly  at  the  limit  turned  out  to  be  stable. 
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13-   Derivation  of  the  sufficiency  of  the  stability 
condition  by  the  energy  method.   We  use  the  notation 
A^u  =  u   .-u.,   A_u  .  =  u  .-u  ._-,  ,   and,  when  applying 
the  summatlon-by-parts  formulae,  we  shall  neglect  terms 
arising  at  the  boundaries.   Then  the  first  equation  of 
(9 '2)  can  be  written 

n+1    n      r  A   n      ,   .,  v  .   n     -, 
u.    -  u.   =   v[A_^Wj_  ^/2-(>'-l)Vj-  1/2^  • 

If  we  always  remember  that  w  and   e  are  centered  at 
j-  1/2  while  u   is  centered  at   J   we  can  drop  the 
suffixes  j,    j  -  1/2  In  the  above.   We  then  multiply 

both  sides  by  u   +u   and  sum  over  (the  suppressed) 

v-|  p       1 1  n  lip 
J.   Denoting  ^Z  (''^•)    ^y     11''^  11   etc.,  we  get 
3  "^ 

II    n+1  II 2        II    ni|2  ,    n+1,    n    .      n,,       ,      -,  w    n+1,    n    .     n. 

u  -      u  =      v(u        +u    ,A,w    )-v(7-l) (u        +u    ,A,e  ) 

+  -r 

/  ^  /  n+1 ,  n,   n^ 
=  -v(A_(u   +u  ),w  ) 

f        -I  \    r  ^    f    n+1 ,  nv  n^ 
+v(7  -1)  (A_(u   +u  ),e  )  , 

after  summing  by  parts.   The  other  equations  of  (9.2) 
give  directly 

n+1 1 2   II  ni|2      ,.  ,n+l  ,n,,n+l. 


w 


-1 1 2   II  n  1)2      ,.   n+1   n,  n+1, 
11   -  II  w  II   =   v(A_u   ,w  +w    )  , 


(7-l)(||e"-^^|2-||e-|h=   (7-l)^(A_A^e-+\e-+e-+l) 

,   ^ ,   , .   n+1   n ,  n+1 , 
-  (7-I)  v(A_u    , e  +e    ) 

By  adding  these  three  equations  together  we  obtain 
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caheellatlonof  most  cross-terms  involving  quantities 
evaluated  at  times   nAt   and  (n+l)At.   We  put 

(15.1)      S^     =     l|u^||2  +||w^||2+   (7-l)||e^l|2  +    ■ 

1/      n\       Ha      ni|2,       ,.      n,      -,\n     n^ 
+  -g(7-l)    |i,iA_|_e    II     +  v(A_u    ,(7-l)e   -w    ); 


then  we  find  that   S^   is  non- Increasing: 


S^.,-S.  =   -  I  (7-1)  ^i[llA^e^||2  +  2(A^e^,A^e^+l) 


'n+1  "n       2  v/  ^.  K-|ii-+-  II   .  ^^"+-  ,^^- 

n+3 


+  II  A  e 


n 


-  ^  0   . 

To  deduce  stability  properties  we  need  to  show  that  S 
Is  a  norm  equivalent  to  the  Lp-norm.   We  have  first  to 
dominate  the  Inner  products  In  S   by  the  strictly 
positive  terms .   It  Is  clear  that 

I  (A  u^,w^)  I   <  a||u^|p  +  -  llw^lP     for  any  a>0  . 

o .  •  T   -1     I  /  A   n  nr  I     o  II  ni|2  ,  1  m  n||2   „   ^  ^ 
Similarly,   I  (A_u  ,e  )|  _<  Pn||u  ||  +-g-||e  ||   ,  p^  >  0. 

But  we  also  have 
|(A_u^e'^)|   =   I(u^,A^e^)|   <  ^  p^  Axl|u'^||2+ ^^llA^e"^  ||2, 

We  combine  these  last  two  estimates  In  the  ratio 
4  :  1-^   where  0  _<  |  _<  1,   and  so  find  that 

(15.2)  S^  >  aj^i^lp  +  a^lw'^li^  +  a  |e^|p  +  a^||A^e^||2 


where  all  the  a.  >  0   If 
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v[a+|(7-l)P^  +1  (1-1)  (7-1)^2  ^^]  ^     1  '   v/a  <  1  , 

(ZzlM  <   ,,.1),    ana  ii^U^     <|(.-1).  . 

Prom  the  last  three,  we  have  a  >  v  ,   ^i  ^   v^  ,   and 
Pp  >  c/cr  (1-| ) ;   substituting  these  Into  the  first  gives 
the  requirement 

V^[l  +  (7-1)^  (e^  +1  (l-O^^)]   <   1 

for  some  ^   for  which  0  ^  ^  jf  1  •   The  square  bracket  is 
maximized  for  ^  =  l/l+2ij,,   which  gives  the  result 

v^(7+2ix)   <   (1  +  2[i)    , 

which  is  the  same  as  (11.2). 

If  this  condition  is  satisfied,   S    forms  a  norm 

n 

and  therefore  the  difference  scheme  is  stable  with  respect 
to  the  norm.   Our  final  step  is  to  show  that  this  implies 
stability  with  respect  to  the  Lp-norm 

T   =  \\^\      +  llw  II  +  ||e  II   ,   which  is  clearly  _<  const.  S 
by  (13.2).   But   S^  <   S^, 


and 


S^  <  (l+7v)|lu^||^  +  (1+v)  llw^lP  +  (7-l)(l+v)  lle^lP 

To  bound  the  last  term,  we  need  to  appeal  to  the  differ- 
ence equation  again 
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|e^  -  n  A_A^  e^ll^  =^  \\e^    -  v  A_  u^H^   , 


1  .e 


<   2  e    +  ov     u 


Hence  S-,  <  const.  (Tq+T^),   and  stability  Is  proved 
for  the  Lp-norm. 
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